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Abstract. The excitation of internal gravity waves by an entropy bubble oscillating in an isothermal atmosphere is 
investigated using direct two-dimensional numerical simulations. The oscillation field is measured by a projection 
of the simulated velocity field onto the anelastic solutions of the linear eigenvalue problem for the perturbations. 
This facilitates a quantitative study of both the spectrum and the amplitudes of excited jr-modes. 
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1. Introduction 

The problem of the excitation of internal gravity waves 
(hereafter IGWs) in stably stratified media is often studied 
in connection with the possible detection of solar g-modes 
(see e.g. the latest attempt of g-mode detection in the 
GOLF data by Gabriel et al. 2002). As an example, two- 
and three-dimensional numerical simulations of penetra- 
tive convection have shown that it is possible to excite such 
waves below the convection zone from the penetration of 
strong downward plumes into the stable radiative zone lo- 
cated below (Hurlburt et al. 1986; Hurlburt et al. 1994; 
Brandenburg et al. 1996; Brumniell et al. 2002; Dintrans 
et al. 2003). 

Furthermore, IGWs play a major role in stellar evo- 
lution as they can transport angular momentum and/or 
chemical elements over long distances in the stellar in- 
terior. This transport mechanism has been proposed to 
explain the quasi-solid rotation profile of the solar core as 
revealed by helioseismology (Kumar et al. 1999) as well as 
the lithium depletion observed in low mass stars (Talon & 
Charbonnel 1998). However, the correct amount of IGWs 
excited by, say, penetrative convection remains still un- 
known as numerical studies were essentially led from a 
qualitative point of view, so that the mode amplitudes 
were not deduced. 

The main goal of this paper is to show that it is possi- 
ble to determine quantitatively the amplitudes of gravity 
waves propagating in a stable zone of a numerical simu- 
lation using the anelastic subspace. This subspace is built 
from the solutions of the associated anelastic linear eigen- 
value problem for the perturbations. The projection of 
the simulated velocity field onto this basis yields the IGW 
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amplitudes. We will present here the application of this 
method to the simple problem of g-mode oscillations in 
an isothermal atmosphere. The dynamics of such an at- 
mosphere is well known (Brandenburg 1988) so that this 
problem is ideal to illustrate and test the validity of our 
method, before applying it to the more difficult problem 
of IGWs excited by penetrative convection, where the at- 
mosphere is in general non-isothermal. 

After presenting the hydrodynamical model describing 
our isothermal atmosphere (Sect.|21l, we apply two classi- 
cal and widely used methods to detect the (7-modes excited 
by an oscillating entropy bubble and show their limita- 
tions (Sect.|2Jl. Next, we introduce our new method based 
on the anelastic subspace and give the analytic forms we 
found for both the eigenfrequencies and eigenvectors of 
the anelastic set of linear equations for the perturbations 
(Sect-EJ. We then apply this formalism to the same simu- 
lation used to test the two classical methods discussed in 
Sect.Oand show that we now have access to both the spec- 
trum and amplitudes of the IGWs (Sect.|SJ). In particular, 
we illustrate the usefulness of our method by studying the 
influence of the box geometry on the mode amplitudes. 
Finally, we conclude in Sect. by introducing the next 
step of this work, which will consist of the application of 
the anelastic subspace method to the detection of IGWs 
that are stochastically excited by penetrating convective 
plumes. 

2. The hydrodynamical model 

2.1. The basic equations 

We consider the two-dimensional propagation of internal 
gravity waves in an isothermal atmosphere consisting of a 
layer of depth d of a perfect gas at temperature To embed- 
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(a): time = 0.0 



ded between two horizontal impenetrable boundaries. The 
fluid is non-rotating and stratified with constant gravity 
and its properties like its kinematic viscosity, y, and spe- 
cific heats, Cp and c„, are assumed to be constant (with 
7 = Cp/c^ = 5/3 for a monatomic gas). Initially, the layer 
is in hydrostatic equilibrium such that its density is given 

by 



(1) 



it points downward in the same 
2 = corresponds to the top 



po{z) = ptopexp(z/iJ), 

where z is depth, i.e. 
direction as gravity g, 
of the box, H — c^/^g is the pressure scale height, 
Cs — VlR*Ta — const denotes the speed of sound (i?* 
being the gas constant) and iV^ = (1 — l/^)g/H = const 
is the square of the Brunt- Vaisala frequency. 

We choose the depth of the layer d as the length unit, 
\fd[g as the time unit and ptop and Tq as the units of 
density and temperature, respectively. The evolution of 
this layer is governed by the equations of conservation of 
mass, momentum and energy 
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where D/Dt ^ d / dt + u ■ W is the total derivative, p the 
density, u the velocity, e the internal energy with P = 
(7 — l)pe, and S is the trace-free rate of strain tensor, 
given by 



S — 



1 



dxi 



V • u 



0.0 
0.2 7 
0.4 
0.6 

0.8 
1 .0 



0.0 

0.4 

0.6 

0.8 
1 .0 



0.0 
0.2 
0.4 
0.6 

0.8 
1 .0 



* 




* 



* 



O.b 



(b): time = 2.5 




0..5 



I 



V 



-O.b 



0.0 



O.b 



Fig. 1. Velocity field superimposed on a grey scale 
representation of the entropy perturbation for a two- 
dimensional simulation of an entropy bubble oscillating in 
an isothermal atmosphere. The asterisks shown in panel 
(a) are used in Fig. 



2.2. Boundary conditions 

When imposing the boundary conditions, we assume the 
top and bottom surfaces to be stress-free with the fixed 
temperature (or, equivalently, internal energy), that is 

dua: du^ 

— — — — — — Uz — V and e = cq on z = U, 1. 

oz oz 

In the horizontal direction, we shall impose periodic con- 
ditions for all fields. 

2.3. Numerics 

The fully nonlinear set of equations l(2Jl is solved using 
the hydrodynamical code described in Nordlund & Stein 
(1990) and Brandenburg et al. (1996). The time advance 
is explicit and uses a third order Hyman scheme (Hyman 
1979). All spatial derivatives are calculated using sixth 
order compact finite differences (Lele 1992). The two- 
dimensional simulations presented here were done with 



the same resolution 256 x 300 (i.e. 256 zones in the hori- 
zontal direction and 300 in the vertical) and aspect ratios, 
A = Lz/Lx, between 2 and 6. Here, and denote the 
domain size in the horizontal and vertical directions, re- 
spectively. Moreover, the dimensionless sound speed and 
the gravitational acceleration are also the same for all 
runs and set equal to Cg = 5 = 1. As a consequence, the 
dimensionless pressure scale height is H = 3/5 (i.e. the 
box extends over 1.66H in the vertical direction so that 
the density contrast between bottom and top is around 
Pbot/ptop — 5.3), while the Brunt- Vaisala frequency is 
N = — 0.82. Finally, the dimensionless kinematic 

viscosity v has been set equal \o v = 10""^ in all simula- 
tions, except in the study of the influence of the viscosity 
on the mode damping rates fSect. l^^ where we also con- 
sidered smaller values down to = 5 x 10~^. 

As initial condition to excite IGWs in the layer, 
we choose the perturbation to be an entropy bubble 
(Brandenburg 1988). For a given entropy perturbation, 
it is indeed easy to find the approximate position of the 
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point around which the bubble will begin to oscillate and 
thus generate IGWs. In dimensionless units we have 



vertical velocities 



0.03 r 



N'^ ^ — — 

dz 



As' 



so Az 



As' 3 , 



where s' = s/cp — (l/7)lnP — In p is the dimensionless 
entropy. All simulations start with an initial (negative) 
entropy perturbation As' = —0.2 located at z = 0.2 so 
that Az = 0.3. This cold bubble then falls down to the 
middle of the layer, z = 0.5, and begins to oscillate around 
this point, generating IGWs. We note that the entropy 
perturbation is done at constant pressure AP = 0, i.e. 
Alnp = - As'. 

Figure ^ shows an example of such a two-dimensional 
hydrodynamical simulation where the velocity field has 
been superimposed onto a grey scale representation of the 
entropy. At t = 0, the negative entropy perturbation looks 
like a bubble at x = and z = 0.2 (Fig.^) and no veloc- 
ity field is present. Under the effect of gravity, the bubble 
begins to descend (Fig. ^) and stabilizes at the depth 
z = 0.5. The bubble thus oscillates around this equilib- 
rium position and IGWs are generated in the whole do- 
main (Fig. nj;). The question is what is the spectrum and 
amplitude of the internal waves excited by this oscillating 
bubble. 



3. Measuring IGWs excited by the oscillating 
bubble from classical methods 

Two techniques are commonly used to measure wave fields 
propagating in direct hydrodynamical simulations: 

Method 1 (hereafter Ml): the simplest method consists of 
firstly, recording the vertical velocity at a fixed point 
and, secondly, performing a Fourier transform of the 
time series. This method was used for example by 
Hurlburt et al. (1986) to detect IGWs in a simulation 
of penetrative convection. 

Method 2 (hereafter M2): a more precise method consists 
of, firstly, taking horizontal Fourier transforms of the 
vertical mass flux pw for every time step and, secondly, 
computing power spectra for the individual time series. 
Peaks corresponding to acoustic or gravity oscillations 
thus appear at a given horizontal wavenumber in 
the (z, w)-plane. This method was used for example 
by Stein & Nordlund (1990) to determine the acoustic 
modes that are excited in their numerical simulations 
of solar granulation. 

In Fig. 121 wc summarize the application of method 
Ml to detect the IGWs propagating in the simulation of 
Fig. n] As expected for wave propagations, vertical ve- 
locities oscillate approximately periodically around zero 
at the three different levels (Fig. |3l). However, when we 
look at the corresponding temporal power spectra, only 
very few modes appear: the fundamental radial acoustic 
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Fig. 2. Detection of IGWs using method Ml. (a): time 
evolution of the vertical velocity recorded at three differ- 
ent positions in the domain (denoted by stars in Fig. QJl) . 
(b): corresponding temporal power spectra. Peaks below 
the buoyancy frequency N — 0.82 correspond to inter- 
nal gravity waves whereas the peak around lu = n is the 
fundamental acoustic mode. 



mode^ and one or two gravity modes below the bounding 
Brunt- Vaisala frequency N « 0.82 (Fig. We recall 
that gravity waves cannot propagate along the vertical 
direction so that the g-modes necessarily correspond to 
nonradial modes (e.g. Unno et al. 1989). Since method 
Ml does not take into account the horizontal dependence 
of the modes, nonradial y-modes with almost the same 
frequency but different horizontal wavenumber kx cannot 
be distinguished. This degeneracy in the mode degrees ex- 
plains why only very few g-modes are captured with this 
method. Here and in the following we present the horizon- 
tal wavenumbers in terms of their smallest possible finite 
value, so that kx = (2Tr/Lx)£ where £ denotes the mode 
degree and is an integer, £ = [0, 1, 2, ..]. 

In Fig. 121 we summarize the application of the second 
method M2 to detect IGWs in the same simulation of 
Fig. 121 We now have access to more informations than with 
method Ml. At first, the £ ~ diagram confirms: (i) that 
the mode around a; = tt is indeed a radial one and that 
it corresponds to the fundamental acoustic mode as no 
radial node is visible; (ii) that gravity modes are strictly 
nonradial modes as no peaks are visible below the dotted 
line u = N. 

Second, gravity modes with similar frequency but dif- 
ferent degrees £ are now well separated. As an example, 
the fundamental g-mode at £ ^ 1 and the first overtone at 
£ = 2 have almost the same frequency, uj ~ 0.6. Using the 
(z, [x))-plane for different ^!'s allows us to separate modes 
and to show that one mode corresponds to a fundamental 



^ It corresponds to sound waves making successive motions 
back and forth in the vertical direction with a period T simply 
given by r = 2/cs = 2, so the pulsation is ui = 2n/T = n. 
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Fig. 3. Detection of IGWs using method M2. Left: tem- 
poral power spectra in the (z,w)-plane for three different 
degrees £ = [0,1,2]. Right: the resulting spectra after an 
integration in depth. Dotted lines shown in all plots cor- 
respond to the bounding frequency uj = N ^ 0.82, above 
which g-modes cannot exist. 

one while the other one is a first overtone [note the two 
"bumps" around to ~ 0.6 in the (z, a;)-diagram for £ — 2]. 

4. A new method using the anelastic subspace 

The classical methods Ml and M2 presented above have 
some disadvantages which render them inappropriate for a 
careful detection of IGWs in hydrodynamical simulations: 

— the main disadvantage of method Ml is its k^- 
degeneracy, that is, the horizontal dependence of the 
modes is not taken into account. As a consequence, we 
have only access to a rough spectrum where only very 
few peaks appear. 

— method M2 takes into account this horizontal depen- 
dence so that more modes are detected. But what 
about modes amplitudes? For example, this method 
does not permit to quantify the amount of kinetic en- 



ergy due to IGWs propagating in the box so that only 
qualitative information, such as the knowledge of the 
spectrum of excited g-modes, is possible. 

Finally, we recall that the goal of this work is to find a 
tool allowing us to measure quantitatively IGWs that are 
stochastically excited by penetrative convection. In view 
of this it is clear that methods Ml and M2 are not well 
adapted to this problem as Fourier transforms done to find 
the spectrum are calculated over the whole simulation. 

Our new method, which eliminates these shortcomings, 
is inspired by the work of Bogdan, Cattaneo & Malagoli 
(1993, hereafter BCM) who have developed a tool to mea- 
sure the acoustic emission generated in their simulations of 
turbulent convection. They extracted the acoustic field by 
projecting the convection field onto the acoustic subspace 
build from the eigenfunctions of the associated linear oscil- 
lations problem. Here we are interested in gravity waves 
rather than sound waves and this allows us to adopt a 
simpler procedure than theirs. Thus, we project our sim- 
ulation data onto the anelastic subspace, that is, we filter 
out the acoustic waves in the linear system of oscillation 
equations. For a time dependence of normal modes of the 
form expiiut), and in the ideal limit 1^ = 0, the anelastic 
set of equations reduces to (Dintrans & Rieutord 2001; 
Ricutord & Dintrans 2002) 



< div(po^) = 0, 
J, = for z = 0,l. 



(3) 



where ^ = {^x,^z) denotes the Lagrangian displacement 
(the associated velocity field being v = iuj^), U = P' /po 
the eulerian perturbation in reduced pressure and po{z) 
the equilibrium density profile. 

Because we adopt periodic boundary conditions in the 
x-direction, we can seek solutions of the form 

^xix,z) = £,xiz)cos{k^x) and ^^{x, z) = ^^{z) sin{kxx), 

such that the anelastic system (jSJ reads after the elimina- 
tion of n 



= 



dz 

= for 



dz 
z = 0, 1. 



(4) 



This last system can be formally rewritten as a generalized 
eigenvalue problem of the form 



(5) 



where ipin = {£,x,S,z)'^ is the eigenvector associated with 
eigenvalue and Ma, -^b denote two differential op- 
erators. 
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As shown in Appendix^ the eigenfunctions xjjgn which 
are solutions of the eigenvalue problem Q are orthogonal, 
that is we have 

f-i 



(6) 



for normalized eigenfunctions (here the symbol f denotes 
the Hermitian conjugate). As a consequence, the set i/j^n 
forms an orthogonal basis onto which the simulated veloc- 
ity field V can be projected, i.e. 

+CX; +00 COs{kxX) 
v{x,Z,t) =^^{4len,v)lpen{z) +Vrcst, (7) 

1=1 n=0 sm{kxx) 

where the i-sum begins at £ = 1 as no gravity modes prop- 
agate radially, i.e. — {2tt/Lx)£ should be non-zero (see 
Sect.OJ. In this decomposition, the first term corresponds 
to the IGW contribution while all other rates, such as the 
bubble displacement or the acoustic waves, were collected 
in a term which we identify hereafter as "rest" . Finally, 
following BCM, we identify the amplitude of each gravity 
mode by the time-dependent coefficient 

C£n{t) = {4>£n,v). (8) 

In the case of an isothermal atmosphere, for which po 
and N are simply given by analytic expressions for 
both eigenvalues and eigenvectors of the anelastic system 
(0J can be found (see Appendix 0) and we give here only 
the result: 

-1-1/2 

1 + ^1 + kt 



1 



ten 



C_ 



exp(-z/2iJ) 



fcz coa{k^z) + — sm{k^z] 



C exp(— z/2iJ) sin(fc2z). 



where we have introduced the vertical wavenumber k^ — 
[n + l)7r. Because of these analytic expressions, the con- 
struction of the anelastic basis used in the projection {7)) 
is straightforward. However, in practice, we do not project 
directly onto the two-dimensional velocity field but rather 
its horizontal Fourier transform at a given k^ , so that the 
procedure is the following: 

1. given a fca;-value, we first deduce the anelastic eigenval- 
ues uin and eigenvectors i/'fn = (^^"jCz")"^ from their 
analytic expressions ©. 

2. we build the anelastic subspace at this horizontal 
wavenumber k^ with, say, the first five anelastic modes 
n= [0,4]. 

3. we Fourier transform the simulated velocity field 
v{x^ z, t) to obtain the new complex-type field vi{z, t). 

4. we project vi{z,t) onto the anelastic basis, built at 
step 2, and deduce the mode amplitudes (|SJ for each 
order n. 

We will now illustrate the efficiency of our IGW de- 
tection method by measuring IGWs excited in the same 
simulation used to present methods Ml and M2 in Sect.|21 





depth 



Fig. 4. Vertical structure of the fundamental anelastic 
mode (dark line) and first overtone (grey line) a,t £ = 1 
(a) and £ — 2 (h). Solid lines mark the horizontal dis- 
placement while dot-dashed lines are for the vertical 
one ^z. Eigenfunctions have been normalized by imposing 
Jo \'4^in\^Podz = 1. 



(9) 5. Results 

5.1. The anelastic eigenvectors: comparisons witli tfie 
simulation profiles 



Figure^lshows four anelastic modes at ^ = 1 (a) and £ — 2 
(b) with n — (dark lines) and n = 1 (grey lines). We 
remark that anelastic eigenfunctions hardly change with 
the ^- values as does not explicitly depend on k^, while 
only the ^^j-amplitudc involves k^. As a consequence, it is 
the ratio £,z/£,x oc k^ which makes sense so that S,z becomes 
more and more important compared to with increasing 
kx , as observed in Fig. ^ when comparing solid to dot- 
dashed lines for the same order n. In the case of method 
M2, the use of the vertical mass flux pqVz (instead of, say, 
PqVx) is justified a posteriori as a good proxy for detecting 
IGWs propagating in the simulation (see Fig.OJ. 

Once the anelastic eigenfunctions are known, one can 
compare the corresponding vertical mass flux pov^ with 
the ones deduced from the power spectra used in method 
M2. Indeed, the peaks appearing at a given frequency in 
Fig. |3] correspond to g-modes so that their vertical profiles 
can be related to the eigenfunctions of the linear modes. 
This was done by BCM for the case of acoustic modes 
propagating in their convection simulations. 

The upper panel of Fig. [3 shows the depth dependence 
of the two peaks appearing at £ = 1 in Fig. El which corre- 
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depth 

Fig. 5. Upper: same as Fig. Oat i = I, except that am- 
plitudes have been restored using third direction. Lower: 
comparisons between the normalized vertical mass flux 
poVz profiles deduced from the upper power spectra (solid 
lines) with the ones obtained from the anelastic modes 
shown in Fig. 0]i (dotted lines). 



spond to the n = (wio = 0.568) and n = 1 (wn = 0.363) 
modes, respectively. We thus recover what BCM called the 
"shark fin" profiles, that is, the power at any depth is lo- 
cated around the theoretical anelastic frequency. Taking 
a mean profile around eigenfrequencies loiq and ujn al- 
lows us to compare the vertical mass flux observed in the 
simulation with those calculated from the linear anelastic 
modes. The lower panel of Fig.Oshows that the agreement 
between the two profiles is quite remarkable, meaning that 
the anelastic modes reproduce well the long-period dy- 
namics of the oscillating isothermal atmosphere. 
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Fig. 6. Left: real parts of the complex mode amplitude 
Cin{t) for the three modes shown in Fig. ^ Right: corre- 
sponding temporal power spectra. Dotted lines mark the 
location of the anelastic eigenmodes wio = 0.568, uju = 
0.363 and W21 = 0.575, respectively. Amplitudes on each 
plot have been multiplied by one thousand. 

5.2. Time evolution of the mode amplitudes cen{t) 

The left hand panel of Fig. |H| shows the real parts of the 
complex amplitudes cio, cn and C21, that is, the result 
of the projection of the simulation shown in Fig. ^ onto 
the three anelastic modes plotted in Fig. 0] By taking a 
Fourier transform of these sequences (Fig. El right), one 
obtains three peaks which agree remarkably well with the 
theoretical anelastic frequency given by @, i.e. the com- 
plex amplitudes q„ behave as 

Q„ (X exp{iujent). (10) 

This relation means that when an eigenmode is excited in 
the simulation, its corresponding displacement is periodic 
with a period Tg^ = 27r/a;£„. 

As we have access to the real and imaginary parts of 
Cf„, we now compute the amplitude modulus |q„|; Fig. [T] 
shows its time evolution for the three modes in Fig. El A 
regression analysis of the curves max(|c£„|) = f{t) allows 
us to determine the shape of the envelope and we found 
that \cin\ follows a simple exponential decay law of the 
form exp{—at). Because of the periodic behaviour (|10|l . it 
means that the mode amplitude behaves like 

cen oc exp(-ai) exp(iw£„t), (11) 

where a is a constant which is related to the damping 
process, here the viscosity. Using runs with different kine- 
matic viscosities we have verified that a (x u, i.e. every 
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250 



Fig. 7. Time-evolution of the amplitude modulus for 
the three modes shown in Fig.^ Grey lines superimposed 
correspond to an exponential decay of the form |q„| oc 
exp(— ai). 




Fig. 8. Dependence of the time-averaged kinetic energy 
on degree ^ and radial order n. 



that a Parseval-Bessel relation also exists in the anelastic 
subspace as 



„ +00 

/ p{)v'^{x, z)(\xi\z =y ^ \ctn\^ + rest, 



(12) 



where, as before, the rest term contains all the contribu- 
tions which are not due to IGWs. 

This relation illustrates the usefulness of the anelastic 
subspace when one wants to perform a detailed study of 
the IGW excitation. By using the classical Parseval-Bessel 
relation applied in Fourier space, it is already possible to 
find the amount of kinetic energy contained at a given 
wavenumber kr, since 



p{)v'^{x, z) dx dz 



„ +00 

/en 



(13) 



For example, the dependence of the time-averaged kinetic 
energy on the degree shown in Fig. |H1 demonstrates 
that half of the kinetic energy in the simulation in Fig. ^ 
is contained in the f = 1 degree (the solid line denoted 
by "g- modes + rest"). Unfortunately, this curve gives no 
information about the vertical dependence of this distribu- 
tion so that contributions coming from nonradial acoustic 
modes as well as gravity modes can be present. 

The main advantage of our anelastic subspace method 
is that it lifts this degeneracy by isolating contributions 
that come from g-modes only. Comparing to its classi- 
cal form l|13|l . the anelastic Parseval-Bessel relation H12|l 
introduces the radial order n so that one can extract the 
contribution of the g-mode (£, n) in the kinetic energy and 
then deduce which modes contribute most to it, meaning 
the modes that are the most excited by the oscillating 
bubble. This is what we did in Fig. |H1 where we also plot- 
ted, for each £, the contributions coming from the n = 0, 
n = 1 and n = 2 y-modes. It allows us to show that the 
t = \ contribution is in fact almost entirely composed of 
the (1,0) and (1, 1) g-modes, as they respectively contain 
around 33% and 15% of the time-averaged kinetic energy 
of this simulation. This kind of information is particularly 
relevant when one deals with the problem of wave excita- 
tion. 



mode obeys the same law as that of a linear damped har- 
monic oscillator. Such a relation between the damping of 
the mode and the viscosity is hardly surprising as it can 
easily be deduced from work integrals; see Appendix lO 

5.3. Contributions of g-modes to the time-averaged 
kinetic energy 

We will now show that one of the main advantages of 
the anelastic subspace resides in the fact that it allows to 
calculate the contribution of every g-mode to the time- 
averaged kinetic energy. We demonstrate in Appendix IdI 



5.4. Effects of the box geometry on modes amplitudes 

In order to assess the restriction imposed by the finite size 
of the box, we now study the influence of the box geome- 
try onto the wave excitation. We therefore perform three 
runs with three different aspect ratios (2, 4 and 6) and 
calculate the g-mode contributions to the time-averaged 
kinetic energy in each case. 

Figure El shows the spectral power in the {kx,Lo) plane 
for the three different aspect ratios. One recovers clearly 
the well-known dispersion relation for g-modes (e.g. Stein 
& Leibacher 1974); see also its anelastic counterpart in 
Eq. jnj. One also verifies that increasing the horizontal ex- 
tent of the box allows the power to be distributed among a 
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larger number of horizontal modes. This phenomenon has 
also been observed in the problem of p-modes excitation by 
turbulent convection (Stein & Nordlund 2001; Nordlund 
& Stein 2001). Unfortunately, such diagrams are gener- 
ally not easy to interpret when the system is not isother- 
mal and the modes stochastically driven; see for example 
the three-dimensional simulations of Brandenburg et al. 
(1996) where the resulting {kx,oj) turned out to be very 
noisy. 
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Fig. 9. (fcx — diagrams which three different aspect ra- 
tios (2, 4 and 6). Crosses correspond to the anelastic fre- 
quencies given by 



Contrasting this with the anelastic subspace method, 
it is straightforward to extract at a given wavenumber 
kx the energy from modes with different radial order n. 
As expected, the spectral distribution of mode energies 
integrated over all n collapses onto an universal curve 
(Fig. 110(1 . It turns out that the mode energy scales as 



-2/3 



for small enough values of k^, i.e. for kxLx/A<10. 



However, both the exponent and the cut-off are not unique 
and depend on the size of the initial entropy bubble, as 
verified by running several simulations with bubble radii 
ranging from i? = 0.1 to i? = 0.4. We found that the larger 
the bubble, the smaller is the cutoff value of k^, meaning 
that big bubbles preferentially generate g-modes at large 
scales. It would be interesting to investigate in more detail 



Fig. 10. Distribution of the energy contained in 5-modes 
for different aspect ratios (2,4 and 6) in a log-log repre- 
sentation. 

the origin of such behaviour, but this is clearly beyond the 
scope of this paper. 

6. Conclusions 

We have presented a new and accurate method to mea- 
sure the internal gravity waves propagating in a numerical 
simulation. This is of prime importance when one deals 
with the problem of the generation of IGWs in a stably 
stratified medium such as the radiative zone of a star. The 
knowledge of the vigour of the velocity field in such a zone 
is the main input of any theory of wave transport of angu- 
lar momentum and/or chemicals in stellar radiative zones 
(e.g. Schatzman 1993). 

However, we have shown that classical methods based 
only on Fourier transforms in space and/or time of the 
simulated velocity field do not allow a quantitative deter- 
mination of the mode amplitudes. The reason is mainly 
that the contribution of the propagating g-modes to the 
total energy in the simulation is diluted among the other 
ones due to p-modes and turbulent motions. 

Our method is based on the projection of the simulated 
velocity field onto the anelastic subspace built from the 
solutions of the linear problem for the perturbations. Once 
the resulting time-dependent coefficients of the velocity 
field are computed, one can deduce the amplitudes of every 
g-mode and measure their kinetic energy at every timestep 
of the simulation. 

We have applied this formalism to the measurement 
of IGWs excited by an entropy bubble oscillating in an 
isothermal atmosphere. We have shown that the mode 
amplitudes follow the same law than those of a linear 
damped harmonic oscillator, that is, the periodic ampli- 
tudes decrease exponentially with time, with a time scale 
that is simply proportional to the inverse of the viscosity; 
see Eq. (|11() . By using a Parseval-Bessel relation valid in 
the anelastic subspace, we have deduced the contribution 
of every g-mode to the time-averaged kinetic energy and 
have shown that large-scale motions are mainly composed 
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of n = and n = 1 modes. Such a result is in general 
impossible to obtain in a more complicated settings using 
a classical (fc^, — lu) diagram. 

The fact that the mode amplitudes in the anelastic 
basis follow the same law than those of a linear damped 
oscillator is important in the case of a stochastic excitation 
of IGWs, such as the one encountered when one deals with 
numerical simulations of penetrative convection. Indeed, 
in this case the mode amplitude q„ evolves either chaot- 
ically around zero when the mode is not excited, or in 
a periodic fashion following Eq. Hll|l when the mode is 
excited. By computing time-frequency diagrams of that 
time-dependent amplitude, one can extract the time in- 
tervals during which the mode has been really excited. In 
a forthcoming paper, we will present the application of the 
anelastic subspace method we tested here to the problem 
of IGWs stochastically excited by penetrative convection 
(Dintrans et al. 2004). 
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Appendix B: Anelastic oscillations of an 
isothermal atmosphere 

B.l. Derivation of the anelastic eigen frequencies and 
eigenfunctions 

We start from the anelastic set of equations Q 

r N\ _ _ 1 



kx dz 



dz dz 
= for 2 = 0, 1, 



(B.l) 



where, for an isothermal atmosphere, po{z) and A*'^ are given 
by Eq. Q . By eliminating , we arrive at the following second- 
order differential equation for 



d'g. 1 djz 
dz^ H dz 



fc' f 1 - ^ I = 0, 



(B.2) 



Uo) = e.(i) = 0. 



All coefficients of this differential equation being constants, we 
can seek solutions of exponential type oc exp(rz) to obtain 
the characteristic equation 



Appendix A: Orthogonality of the anelastic 
eigenfunctions 

We start from the complete set of equations written for two 
eigenmodes (ti;i,^i) and (0)2,^2) (Dintrans & Rieutord 2001) 

(1), 



div(po4i) = div(/9o|2) = 0, 



(A.l) 



(1) 



(2) 



for 2 = 0, 1, 



where 11 = P' / Po denotes the eulerian fluctuation of the re- 
duced pressure. Multiplying the ^i-equation by {2 and the {2- 
equation by ^1, and subtracting the resulting two equations 
leads to 



{uji — uj2)poii ■ $2 = div(nip()|2) - div(n2po|i). 



(A.2) 



where we have used the relation po^2 ■ VIIi = div(nip()52) — 
Hi div(po^2) = div(nipo42), as the anelastic approximation 
implies div(po^2) = 0. 

Finally, we integrate this last equation over the volume and 
use the Green-Ostrogradsky theorem to obtain 

iujf-ij^2) I Poii-i2AV= [ Uipoi2-dS- [ n2PoCi-dS'.(A.3) 
Jv Js Js 

The applying of periodic horizontal boundary conditions 
and the condition = at z = and 2 = 1 leads to the 
vanishing of the two right hand side terms. The uniqueness of 
the solutions ensures that uji 7^ lu2 so that we conclude that 
the eigenvectors are orthogonal to each other, i.e. 



poll ■ i2dV = 0. 



(A.4) 



H 



= 0, 



whose general solution is of the form 



-4fcf 



(B.3) 



A^ \ 
1-^ ,(B.4) 



leading to the following solution for 
^^(2) = exp(-2/2Ji') Acos 



B sin 



/A 



.(B.5) 



To find the quantization rule for A, we consider the boundary 
conditions on ^z, 




(B.6) 



= (n-|-l)7r or A„ = 4(n -|- 1) tt 



where n = [0,1,2...] is the radial order of the mode. [The 
presence of the n-\-l factor in the quantization rule is consistent 
with the fact that A is strictly non-zero.] We recall that the 
definition IB. 41 for A allows us to deduce an expression for the 
anelastic eigenfrequencies. 



OJln — A 



1 + ^(4^ + '^^ 



-1/2 



with kz = {n + 1)-K, (B.7) 



and their associated eigenfunctions 



= f exp(-2/2^) 



kz cos{kzz) + — — sin{kzz) 



(B.i 



^i" = Cexp(-2/2H) sm{kzz). 



The remarkable fact to be noted here is that the vertical la- 
grangian displacement ^f" does not depend on kx (the horizon- 
tal wavenumber) but only on the mode order n (see Fig. |1J. 
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Table B.l. Exact eigenfrequencies w for the first five g- 
modes of degree £ — 1, corresponding to an horizontal 
wavenumber = 217/ = tt (here L^. = 2), in the com- 
plete case, in the anelastic approximation, and the corre- 
sponding relative errors. 



Order n 


UJ 




rel. error 





0.572 054 


0.567 455 


0.804% 


1 


0.363 084 


0.362 606 


0.132% 


2 


0.257 381 


0.257 295 


0.033% 


3 


0.197 644 


0.197 621 


0.012% 


4 


0.159 920 


0.159 912 


0.005% 




.0 I 

0.0 



dep:h 



Fig. B.l. Comparison of the complete (solid lines) and 
anelastic (dashed lines) normalized eigenfunctions for 
€ = 1 and n = [0,1,2]. 



B.2. Comparisons with the complete case 

The set of equations for the complete case is obtained by, first, 
linearizing the system for infinitesimal adiabatic perturba- 
tions and, second, by assuming the time and horizontal depen- 
dences of the modes to be exp(ia;t) and cos{kxx) or sm(k^x), 
respectively. We thus arrive at 



n 



(B.9) 



az \ 9 



L = for z^O, 1. 



Here 11 = P'/Po and p' are the eulerian perturbations in re- 
duced pressure and density, respectively. This system has been 
obtained under the classical Cowling approximation which as- 
sumes that perturbations in the gravitational potential are neg- 
ligible (Cowling 1941). We refer the reader to the book of Unno 
et al. (1989) for a complete derivation of this system in the case 
of a spherical geometry. 

Lamb (1909, 1932) found analytical solutions of the system 
JB.9II in the case of an isothermal atmosphere. By eliminating 



^a; and n in the system IIB.9t , we obtain a second-order differ- 
ential equation for only. 



dz^ H dz 



I ?.(0) = C.(1) = 0. 



.2 ^ ,,,2 



(B.IO) 



Comparison with Eq. ljB.2^ shows that the anelastic approxi- 
mation applies provided that: 

— the Lamb frequency term /(?s is negligible, implying that 



(B.ll) 



which is a condition naturally fulfilled by high-order g- 

modes with ui <^ N . 

the ratio kxCs/N is very large^ 



The solutions of the complete system IIB.lOt can be de- 
duced from the anelastic ones as only the A-term is changing, 
that is, 



^compl — ^ 



(B.12) 



while the same quantization rule applies due to the identical 
boundary conditions Cz(0) = ^z(l) = 0, that is. 



Acompl = 4(n + 1)V^ = Akl. 



(B.13) 



We thus obtain the following quadratic dispersion relation for 
the complete eigenfrequencies 



(^kl + kl + ^-i^ j Lj^ + klclN^ = 0, 



(B.14) 



which is the same than that in Stein & Leibacher (1974). In 
the same way, the eigenvectors of the complete case can be 
deduced from the anelastic ones and we find 



C^ = f exp(^^/2//) 



fcz cos{kzz) + T7 ( - - 7: I sm{kzz) 
H \"f 2 ' 



Cf" = Cexp{~z/2H)sm{k,z) 



Despite its different A-term (|BTT2J, we note that in the com- 
plete case the vertical displacement ^z is the same as in the 
anelastic one, which is simply a consequence of applying the 
same quantization rule IIB.1311 . 

We summarize in Table IB.H the eigenfrequencies for the 
complete case deduced from IIB.14t for the first five g-modes 
of degree ^ = 1 and compare them with their anelastic coun- 
terparts given by the formula IIB.711 . The agreement is quite 
remarkable since the relative error done by the anelastic ap- 
proximation is less than 1% from the fundamental g-mode. 
Concerning the associated eigenvectors, we showed that the 
vertical displacements ^z are the same in both cases so that 
only the horizontal ones differ. However, Figure lHrl shows that 
the agreement on the normalized ^^-eigenfunctions is also very 

^ This last condition implies a clean separation between 
the acoustic and gravity spectra such that a mixed "gravito- 
acoustic" mode is not possible, i.e. the two spectra do not over- 
lap (see Rieutord & Dintrans 2002). 
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except for tlie mode at n = for wliicli tlie condition 
||rTT)i is not well fulfilled. 

In conclusion, the filtering of acoustic waves from the dy- 
namics of the isothermal atmosphere does almost not change 
the p-mode eigenfrequencies uie„ and eigenfunctions tpin ~ 
i^x,^z)'^- The basis built from these anelastic eigenfunctions 
can therefore be used with good confidence to project out ve- 
locity fields from hydrodynamical simulations. 

It should be pointed out that the application of the anelas- 
tic subspace method will only give good results when the 
acoustic and gravity spectra are well separated, that is, when 
the ratio k^Cs/N is large. This is the case for the chosen 
isothermal setup, where the resulting anelastic eigenfrequen- 
cies/eigenvectors are very close to those of the complete prob- 
lem (see Table IHtI and Fig. IB.lt . However, this ratio decreases 
with larger stratification, in which case one may be forced to 
solve numerically the oscillation equations lB.9t of the com- 
plete problem instead of those of the anelastic subset (IB.IL 
which does not however introduce any particular difficulty. 



Appendix C: Damping rate of the modes using 
work integrals 

We show in this appendix that it is possible to find the rela- 
tion governing the damping rate of a mode and the viscosity 
using the work integrals formalism. We start from the anelastic 
equations written for the velocity 



' All = - Vn - N'^izez + v^u, 
div(pou) = 0, 



(C.l) 



where \ = a iuj has a non-zero real part due to the vis- 
cous dissipation. We note that the viscous term is normally 
more complicated as there may be a vertical variation of the 
dynamical viscosity jj. = pv [see the general form of this term 
in Eq. (|5J]. However, we take here a simpler form as our aim 
is just to formally prove that the damping rate a is linearly 
related to the viscosity v. 

We multiply the momentum equation by pou* to obtain 



A similar relation also applies in the complete case, pro- 
vided one includes the contribution of pressure fiuctuations to 
the energy, that is. 



Pou* ■ Aitdy 



viscous dissipation 
total energy ' 



where the total energy contains the kinetic {pQV?), potential 
(N'^po^z) and internal (poTl/Cg) contributions. Such relations 
are useful from a numerical point of view as they allow to check 
the accuracy of the computed eigenvalues from their associated 
eigenvectors. 

Appendix D: Parseval-Bessel relation in the 
anelastic subspace 

In this appendix, we will show that it is possible to relate the 
total kinetic energy embedded in the box to the one calculated 
in the anelastic subspace. To do that, we need a similar relation 
than the Parseval-Bessel one which states that the total energy 
is conserved when one works in Fourier space, i.e. 



f{x,y,z)dV^ I 



(D.l) 



where / denotes the Fourier transform of /. 

In our case, we only perform a 1-D Fourier transform of 
the velocity field in the horizontal direction (corresponding to 
a wavenumber kx — {2n/Lx)(.), so that the Parseval-Bessel 
relation for the total kinetic energy at a given time is 



/ pov^{x,z)dxdz ^ / \ve\^podz, 
Jv J z , 



(D.2) 



where vg denotes the velocity field at a wavenumber k^- We 
now introduce the anelastic subspace by using the relation 
which links the velocity field in Fourier space to the one in 
the anelastic space as 



vi{z) = ''^{tpen,Vi)tpeniz) + rest. 



(D.3) 



Xpou^ = — V(pHu') — N^po^zul + vpQU* ■ Ai 



(C.2) 



where we have used dw{pou*) = 0. As in Appendix El we 
finally integrate this last equation over the volume and use the 
Green-Ostrogradsky theorem to obtain 



A f {pou'' + N^po^z) dV = u f pou* ■ AudV, 



(C.3) 



where we also used it* = A^*. As we only have real integrals 
on both sides, we can deduce the following relation for the real 
part a 



pou* ■ AudV 
' J^ipou^ +N^Poez)dV 



(C.4) 



The z-integral in Eq. 1D.2II can then be calculated in the anelas- 
tic subspace as 



\ve\ po dz = 



v\vipo dz. 



= '^^{ipiii.,vi) I v\^pln(z)podz + xest, 

n ^ 

= ^^{if'in,ve){ve,if'tn) -I- rest. 



(D.4) 



+ 00 

\ctn\ +rest. 



The damping rate of a mode is therefore linearly related to the 
viscosity, as found numerically in lllll . 



This last equation means that the total kinetic energy due to 
(7-modes at a given kx and time t is simply given by the sum 
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over the radial orders n of the squares of the mode amplitudes. 
By taking into account the fc-integral, we finally obtain the 
Parseval-Bessel relation written in the anelastic subspace as 



X 



pov^{x, z) AxAz = \cin\^ + rest. (D.5) 
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